Food Insecurity in the United States
Nawar Naseer
Overview
Food insecurity, or a “lack of consistent access to enough food for an active, healthy life”, is a major public health issue in the United States, where approximately 11% of households experienced food insecurity in 2018 (1). The United States Department of Agriculture (USDA) publishes a Food Environment Atlas, which contains data on over 275 variables with respect to food insecurity, including access to food retail options, availability and use of government food assistance programs, health and socioeconomic indicators, etc. These data are collected from various sources, span several years and geographic levels, and have been made publicly available. The goal of this project is to map food insecurity in the United States. Specifically, this report demonstrates:
- how the food retail environment varies across the United States
- how different factors such as race, income, age, food assistance benefits, and health outcomes change with respect to food insecurity across the United States
Geospatial differences are analyzed, with a specific focus on Philadelphia in comparison to
- other counties in Pennsylvania, and
- Baltimore, MD, another large urban county recommended as a suitable comparison by Moriah Hall, a faculty consultant.
Introduction
In the United States, approximately 11% of households were food insecure in 2018 (1). In Philadelphia, this rate is nearly double the national average at 18.3%, and more than one out of every five Philadelphian struggles with basic food access (2). Food insecurity contributes to the development of chronic diseases such as heart disease and diabetes, both of which are two of the United State’s top ten leading causes of death (3). Additionally, food insecurity affects children’s cognitive development, and their ability to excel in school (4). Thus, food insecurity has long-lasting impacts on the productivity of a population, and understanding the regional factors that contribute to food insecurity and adverse health outcomes is vital for creating effective local and national interventions.
A food secure environment provides both physical and financial access to nutritious food. By analyzing the food retail environment in context of socioeconmic factors and government assistance programs, this project spans the disciplines of sociology, economics, geography, public policy, and public health. Based on consulting with Dr. Allison Karpyn, an analysis on a vulnerable but previously under-researched population, senior citizens, has also been incorporated by assessing the impact of age on food access. Another consultant, Moriah Hall, suggested doing a comparison of Philadelphia, PA to Baltimore, MD, a peer urban county. This project thus aims to use a multidisciplinary approach to understanding region-specific factors influencing food insecurity, and in doing so, may provide critical insight on how to tackle food insecurity in the United States.
Methods
Downloading data and necessary packages
The main data set used for this project was obtained from the USDA’s publicly available data on the Food Environment Atlas. The current version, last updated March 27, 2018, was downloaded directly from the website. The downloaded Excel file contained several sheets pertinent to this project and have been loaded into R as a list of dataframes called “mysheets”. These dataframes contain county-specific information on access and proximity to food retail providers (grocery stores, convenience stores, farmers’ markets, etc), food assistance usage, food insecurity rates, diabetes and obesity rates, and socioeconomic indicators including race, age, income, etc.
To create maps, the counties polygon dataset from BMIN503 Assignment 5 has been downloaded and used.
All necessary packages will be installed and loaded first.
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.2.1 ✔ purrr 0.3.3
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(dplyr) #used to clean data
library(ggplot2) #used to generate plots, static maps
library(RColorBrewer) #used to color plots, maps
library(sf) #used for leaflet to work properly## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0The data downloaded has been saved into the working directory and is loaded into R. Note that names of some of the sheets were slightly modified in Excel to exclude spaces and make them R-compatible.
#Load in data
#The following function allows one to read in multiple sheets from one Excel file into a list, with each sheet representing a separate dataframe within the list.
read_excel_allsheets <- function(filename, tibble = TRUE) {
sheets <- readxl::excel_sheets(filename)
x <- lapply(sheets, function(X) readxl::read_excel(filename, sheet = X))
if(!tibble) x <- lapply(x, as.data.frame)
names(x) <- sheets
x
}
mysheets <- read_excel_allsheets("USDAFoodEnvironmentAtlasMarch2018.xls")
head(mysheets$ACCESS) #to check that data has been loaded properlyLoad in the counties polygon dataset to use for creating maps.
The following themes will be used throughout the report for generating maps.
#For static maps:
myPalette <- colorRampPalette(brewer.pal(9, "YlOrRd"))
my_theme <- function() {
theme_minimal() +
theme(axis.line = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
panel.grid = element_line(color = "white"),
legend.key.size = unit(0.8, "cm"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
plot.title = element_text(size = 12))
}Cleaning data
The next step is to clean the data using the dplyr and tidyverse packages. Each dataframe will be updated to only include variables that are of interest to us for this project.
Access <- mysheets$ACCESS %>%
select("FIPS", "State", "County", "PCH_LACCESS_POP_10_15",
"PCT_LACCESS_POP10",
"PCT_LACCESS_POP15", "PCH_LACCESS_LOWI_10_15",
"PCT_LACCESS_LOWI10",
"PCT_LACCESS_LOWI15", "PCH_LACCESS_HHNV_10_15",
"PCT_LACCESS_HHNV10",
"PCT_LACCESS_HHNV15", "PCT_LACCESS_SNAP15",
"PCT_LACCESS_CHILD10",
"PCT_LACCESS_CHILD15",
PCH_LACCESS_SENIORS_10_15,
PCT_LACCESS_SENIORS10,
PCT_LACCESS_SENIORS15,
"PCT_LACCESS_WHITE15",
"PCT_LACCESS_BLACK15",
"PCT_LACCESS_HISP15",
"PCT_LACCESS_NHASIAN15",
"PCT_LACCESS_NHNA15",
"PCT_LACCESS_NHPI15",
"PCT_LACCESS_MULTIR15")
Stores <- mysheets$STORES %>% #Turns out I don't need to add "" for the column names!
select(FIPS, State, County, GROCPTH09,
GROCPTH14,
PCH_GROCPTH_09_14,
SUPERCPTH09,
SUPERCPTH14,
PCH_SUPERCPTH_09_14,
CONVSPTH09,
CONVSPTH14,
PCH_CONVSPTH_09_14,
SPECSPTH09,
SPECSPTH14,
PCH_SPECSPTH_09_14,
SNAPSPTH12,
SNAPSPTH16,
PCH_SNAPSPTH_12_16,
WICSPTH08,
WICSPTH12,
PCH_WICSPTH_08_12)
Restaurants <- mysheets$RESTAURANTS %>%
select(FIPS, State, County,
FFRPTH09,
FFRPTH14,
PCH_FFRPTH_09_14,
FSRPTH09,
FSRPTH14,
PCH_FSRPTH_09_14)
Assistance <- mysheets$ASSISTANCE %>%
select(FIPS, State, County,
REDEMP_SNAPS12,
REDEMP_SNAPS16,
PCH_REDEMP_SNAPS_12_16,
PCT_SNAP12,
PCT_SNAP16,
PCH_SNAP_12_16,
SNAP_PART_RATE08,
SNAP_PART_RATE13,
PCT_NSLP09,
PCT_NSLP15,
PCH_NSLP_09_15,
PCT_FREE_LUNCH09,
PCT_FREE_LUNCH14,
PCT_REDUCED_LUNCH09,
PCT_REDUCED_LUNCH14,
PCT_SBP09,
PCT_SBP15,
PCH_SBP_09_15,
PCT_SFSP09,
PCT_SFSP15,
PCH_SFSP_09_15,
REDEMP_WICS08,
REDEMP_WICS12,
PCH_REDEMP_WICS_08_12,
PCT_WIC09,
PCT_WIC15,
PCH_WIC_09_15,
PCT_CACFP09,
PCT_CACFP15,
PCH_CACFP_09_15)
Insecurity <- mysheets$INSECURITY
Price <- mysheets$PRICES_TAXES %>%
select(FIPS, State, County, MILK_PRICE10)
Local <- mysheets$LOCAL %>%
select(FIPS, State, County,
FMRKTPTH09,
FMRKTPTH16,
PCH_FMRKTPTH_09_16,
PCT_FMRKT_SNAP16,
PCT_FMRKT_WIC16,
PCT_FMRKT_WICCASH16,
PCT_FMRKT_SFMNP16)
Health <- mysheets$HEALTH %>%
select(FIPS, State, County,
PCT_DIABETES_ADULTS08,
PCT_DIABETES_ADULTS13,
PCT_OBESE_ADULTS08,
PCT_OBESE_ADULTS13)
SE <- mysheets$SOCIOECONOMIC %>%
select(FIPS, State, County,
PCT_NHWHITE10,
PCT_NHBLACK10,
PCT_HISP10,
PCT_NHASIAN10,
PCT_NHNA10,
PCT_NHPI10,
PCT_65OLDER10,
PCT_18YOUNGER10,
MEDHHINC15,
POVRATE15,
PERPOV10,
CHILDPOVRATE15,
PERCHLDPOV10)Further analyses will be conducted in the Results section.
Results
To address food insecurity, interventions need to address three main facets: 1. availability of food: Availability of healthy, nutritious food is impacted by the physical food retail environment of a neighborhood. Are there grocery stores, farmers’ markets, or other avenues through which residents can acquire food? 2. access to food: Access refers to both physical availability of food, but also the financial ability to acquire the food. How expensive are basic food items? Are there stores that accept government assistance benefits that individuals can use to purchase food? 3. utilization of food: Utilization of food refers to individuals actually purchasing and consuming healthy food items. Although this dataset does not have any comsumption data, we can use redemption of government food assistance benefits as a proxy for utilization, as these food assistance programs have restrictions on the types of food individuals can purchase with them.
The next few sections will assess some of these questions and issues.
Lay of the land: state of food insecurity in the United States
The map below shows the three-year household food insecurity average for each state from 2013-2015. What immediately jumps out is that most states in the US have fairly high insecurity rates of above 10%. The situation is particularly dire the southern United States with rates generally being 15% or higher, with a few exceptions. States with particularly high food insecurity include: Oregon, Arizona, Louisiana, Alabama, Missouri, and Kentucky.
#Create new dataframe with variables of interest
Insec1315 <- Insecurity %>%
select(FIPS, State, County, FOODINSEC_13_15)
Insec1315$GEO_ID <- paste0("0500000US", Insec1315$FIPS)
Insec1315map <- inner_join(counties, Insec1315, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Insec1315 <- paste0(Insec1315map$State, # paste0 to append tract name with other relevant text
"<br>Food Insecurity, Three Year Average (2013 - 2015): ", round(Insec1315map$FOODINSEC_13_15, 1))
# Basic leaflet map
library(leaflet)
leaflet(Insec1315map) %>% addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(FOODINSEC_13_15), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Insec1315) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~FOODINSEC_13_15, # variable to be passed to palette function
title = 'Food Insecurity (%), Three Year Average (2013 - 2015)', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()Food insecurity: Philadelphia, PA versus Baltimore, MD
One limitation of the three-year food insecurity rates data is that it is at the state level, and we cannot compare Philadelphia to other counties within or outside of PA. However, we can compare the state of Pennsylvania to Maryland. The bar graph below shows that the three-year average in PA and MD are very similar, making them ideal for drawing comparisons across other variables. While the food insecurity rate decreased in Maryland over time, the rate did not change in Pennsylvania, and it will be interesting to explore other differences between these two states, and specifically between Philadelphia, PA and Baltimore, MD throughout this report.
#Create dataframe of interested variables
Insec10121315_PM <- Insecurity %>%
select(FIPS, State, County, FOODINSEC_10_12, FOODINSEC_13_15) %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(Insecurity_2010_2012 = FOODINSEC_10_12, Insecurity_2013_2015 = FOODINSEC_13_15) %>%
select(-FIPS)
IPM10121315 <- data.frame(State = rep(c("Maryland", "Pennsylvania"), 2),
Year = rep(c("2010-2012", "2013-2015"), each=2),
Percent = c(Insec10121315_PM$Insecurity_2010_2012, Insec10121315_PM$Insecurity_2013_2015))#Create a bar graph to visualize
IPM10121315plot <- ggplot(data = IPM10121315, aes(State, Percent, fill=Year)) + #Load data, specify variables
geom_bar(stat="identity", position="dodge") + #Add a visual layer that is a barplot
labs(title="Food Insecurity, Three Year Average") + labs(x="State", y="Household food insecurity (%) 3 year average")
print(IPM10121315plot)
Availability of food
In this section, we will explore the food retail environment across the United States. We are using food retail providers such as grocery stores as a proxy to assess availability of food. The maps below depict the number of grocery stores, supercenters and club stores, convenience stores, and farmers’ markets at the county level relative to the population.
Grocery stores
For the purposes of this dataset, grocery stores included supermarkets, smaller grocery stores, and deli stores. Most counties across the United States have 0.5 or fewer grocery stores per 1000 individuals. Generally, the northwestern half of the US has more counties with higher numbers of grocery stores. Most of the eastern half of the United States, by contrast, have fewer grocery stores. This may be one of the factors that contribute to the lower food insecurity rates in the northern half of the United States, and higher food insecurity rates in the southeastern half.
#Create new dataframe with variables of interest
Grocery2014 <- Stores %>%
select(FIPS, State, County, GROCPTH14)
Grocery2014$GEO_ID <- paste0("0500000US", Grocery2014$FIPS)
Grocery2014map <- inner_join(counties, Grocery2014, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Grocery2014 <- paste0(Grocery2014map$County, ",", " ", Grocery2014map$State, # paste0 to append tract name with other relevant text
"<br>Grocery stores/1,000 population, 2014: ", round(Grocery2014map$GROCPTH14, 1))
# Basic leaflet map
leaflet(Grocery2014map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(GROCPTH14), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Grocery2014) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~GROCPTH14, # variable to be passed to palette function
title = 'Grocery stores/1,000 population, 2014', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific grocery stores
Looking specifically at Pennsylvania, Philadelphia has a higher number of grocery stores relative to other counties in Pennsylvania. This may be due to how urban Philadelphia is relative to the other counties in PA.
#Create new dataframe with variables of interest
GroceryPA2014 <- Grocery2014 %>%
filter(State=="PA")
GroceryPA2014map <- inner_join(counties, GroceryPA2014, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorGroceryPA2014mapplot <- ggplot() +
geom_sf(data = GroceryPA2014map) +
geom_sf(data = GroceryPA2014map, aes(fill = GROCPTH14), lwd = 0) +
my_theme() +
ggtitle("County-Level Grocery stores in Pennsylvania, 2014") +
scale_fill_gradientn(name = "Grocery stores/1000 population", colours = myPalette(100))
print(GroceryPA2014mapplot)
Supercenters and club stores
Generally, the eastern half of the US appears to have higher numbers of supercenters and club stores relative to the western half.
#Create new dataframe with variables of interest
Super2014 <- Stores %>%
select(FIPS, State, County, SUPERCPTH14)
Super2014$GEO_ID <- paste0("0500000US", Super2014$FIPS)
Super2014map <- inner_join(counties, Super2014, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Super2014 <- paste0(Super2014map$County, ",", " ", Super2014map$State, # paste0 to append tract name with other relevant text
"<br>Supercenters and club stores/1,000 population, 2014: ", round(Super2014map$SUPERCPTH14, 4))
# Basic leaflet map
leaflet(Super2014map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(SUPERCPTH14), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Super2014) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~SUPERCPTH14, # variable to be passed to palette function
title = 'Supercenters and club stores/1,000 population, 2014', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific supercenters/club stores
Relative to other counties in Pennsylvania, there are fewer superstores and club centers in Philadelphia. Again, this may be due to how urban Philadelphia is. Rural counties may be more likely to have the larger spaces at lower prices that are often needed for large warehouse stores such as superstores and club centers.
#Create new dataframe with variables of interest
SuperPA2014 <- Super2014 %>%
filter(State=="PA")
SuperPA2014map <- inner_join(counties, SuperPA2014, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorSuperPA2014mapplot <- ggplot() +
geom_sf(data = SuperPA2014map) +
geom_sf(data = SuperPA2014map, aes(fill = SUPERCPTH14), lwd = 0) +
my_theme() +
ggtitle("County-Level supercenters and club stores in Pennsylvania, 2014") +
scale_fill_gradientn(name = "Supercenters and club stores/1000 population", colours = myPalette(100))
print(SuperPA2014mapplot)
Convenience stores
Aside from a few exceptions, most of the US is relatively uniform in distribution with respect to the number of convenience stores. Just from looking at the color scales of the maps, tt is interesting to note that there are more convenience stores than any other type of retail option in the United States.
#Create new dataframe with variables of interest
Conv2014 <- Stores %>%
select(FIPS, State, County, CONVSPTH14)
Conv2014$GEO_ID <- paste0("0500000US", Conv2014$FIPS)
Conv2014map <- inner_join(counties, Conv2014, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Conv2014 <- paste0(Conv2014map$County, ",", " ", Conv2014map$State, # paste0 to append tract name with other relevant text
"<br>Convenience stores/1,000 population, 2014: ", round(Conv2014map$CONVSPTH14, 4))
# Basic leaflet map
leaflet(Conv2014map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(CONVSPTH14), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Conv2014) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~CONVSPTH14, # variable to be passed to palette function
title = 'Convenience stores/1,000 population, 2014', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific convenience stores
There are relatively fewer convenience stores in southeastern Pennslyvania. I would have expected urban areas like Philadelphia to have a higher density of convenience stores than rural counties, so this finding is a bit surprising to me.
#Create new dataframe with variables of interest
ConvPA2014 <- Conv2014 %>%
filter(State=="PA")
ConvPA2014map <- inner_join(counties, ConvPA2014, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorConvPA2014mapplot <- ggplot() +
geom_sf(data = ConvPA2014map) +
geom_sf(data = ConvPA2014map, aes(fill = CONVSPTH14), lwd = 0) +
my_theme() +
ggtitle("County-Level Convenience stores in Pennsylvania, 2014") +
scale_fill_gradientn(name = "Convenience stores/1000 population", colours = myPalette(100))
print(ConvPA2014mapplot)
Farmers’ markets
The northwestern United States has relatively higher numbers of farmers’ markets. The counties with high numbers of farmers’ markets are located in states with high and moderate food insecurity rates. Farmers’ markets have been used as an intervention to address food insecurity. It would be intersting to explore whether these counties with high numbers of farmers’ markets were the targets of any interventions, or whether farmers’ markets have organically arisen in these locations.
#Create new dataframe with variables of interest
Farmers2016 <- Local %>%
select(FIPS, State, County, FMRKTPTH16)
Farmers2016$GEO_ID <- paste0("0500000US", Farmers2016$FIPS)
Farmers2016map <- inner_join(counties, Farmers2016, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Farmers2016 <- paste0(Farmers2016map$County, ",", " ", Farmers2016map$State, # paste0 to append tract name with other relevant text
"<br>Farmers' market/1,000 population, 2016: ", round(Farmers2016map$FMRKTPTH16, 4))
# Basic leaflet map
leaflet(Farmers2016map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(FMRKTPTH16), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Farmers2016) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~FMRKTPTH16, # variable to be passed to palette function
title = "Farmers' market/1,000 population, 2016", # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific farmers’ markets
In Pennsylvania, the counties immediately surrounding Philadelphia have fewer farmers’ markets. This may be due to the efforts of The Food Trust and other non-profit organizations, as well as the Department of Public Health, in promoting interventions involving farmers’ markets in Philadelphia.
#Create new dataframe with variables of interest
FarmersPA2016 <- Farmers2016 %>%
filter(State=="PA")
FarmersPA2016map <- inner_join(counties, FarmersPA2016, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorFarmersPA2016mapplot <- ggplot() +
geom_sf(data = FarmersPA2016map) +
geom_sf(data = FarmersPA2016map, aes(fill = FMRKTPTH16), lwd = 0) +
my_theme() +
ggtitle("County-Level Farmers' markets in Pennsylvania, 2016") +
scale_fill_gradientn(name = "Farmers' Market/1000 population", colours = myPalette(100))
print(FarmersPA2016mapplot)
Comparison of Food Retail Environment: Philadelphia, PA versus Baltimore, MD
With the exception of grocery stores, Baltimore has a larger presence of stores serving food than Philadelphia relative to the population. This may account for the lower food insecurity rate observed in Baltimore, MD. However, Philadelphia has nearly double the number of grocery stores as Baltimore relative to population. This may highlight the importance of not only providing availability of food, but also phhysical and financial access to food. Philadelphia’s grocery stores may be less physically accessible, or more expensive than Baltimore’s.
(Note: I have not included farmers’ market in this graph because the data are from 2016, whereas the data for the other retail stores are from 2014.)
#Create dataframe of interested variables
Retail_PHLBAL <- Stores %>%
select(FIPS, County, GROCPTH14, SUPERCPTH14, CONVSPTH14, SPECSPTH14) %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(Grocery = GROCPTH14, Supercenter = SUPERCPTH14, Convenience = CONVSPTH14, Specialized = SPECSPTH14) %>%
select(-FIPS)
RPB <- data.frame(county = rep(c("Baltimore", "Philadelphia"), 2),
Type = rep(c("Grocery", "Supercenter", "Convenience", "Specialized"), each=2),
Number = c(Retail_PHLBAL$Grocery, Retail_PHLBAL$Supercenter, Retail_PHLBAL$Convenience, Retail_PHLBAL$Specialized))#Create a bar graph to visualize
RPBplot <- ggplot(data = RPB, aes(x=Type, y=Number, fill = county)) + #Load data, specify variables
geom_bar(stat="identity", position = "dodge") + #Add a visual layer that is a barplot
labs(title="Presence of Stores in Philadelphia and Baltimore, 2014") + labs(x="Store Type", y="Number of stores per 1000 of population")
print(RPBplot)
Access to food
In this section, we will explore physical and financial access to food.
Physical access to food
The dataset contains information about the number of people with low access to food. Low access in this study is defined as living more than 1 mile from a supermarket, supercenter, or large grocery store if in an urban area, or more than 10 miles if in a rural area. It also contains information about factors that might affect low access, such as age and race.
Low Access to Stores
The map below shows the percent popluation with low access to stores in 2015 across the US. The western United States has higher percentage of individuals with low access to stores. What’s intriguing is that areas with low access to stores do not all overlap with areas of high food insecurity. For example, the northwestern part of the US has relatively lower food insecurity rates, yet there are many counties in this region with a high percentage of the population having low access to grocery stores.
One limitation of this Access dataset is that only large grocery stores and supermarkets were considered, and many individuals may use other retail options, such as convenience stores, to meet their regular grocery needs.
#Create new dataframe with variables of interest
LA2015 <- Access %>%
select(FIPS, State, PCT_LACCESS_POP15)
LA2015$GEO_ID <- paste0("0500000US", LA2015$FIPS)
LA2015map <- inner_join(counties, LA2015, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorLA2015mapplot <- ggplot() +
geom_sf(data = LA2015map) +
geom_sf(data = LA2015map, aes(fill = PCT_LACCESS_POP15), lwd = 0) +
my_theme() +
ggtitle("County-Level Percent Population with Low Access to Stores in the US in 2015") +
scale_fill_gradientn(name = "Percent Population", colours = myPalette(100), limit=range(0, 100))
print(LA2015mapplot)
PA-specific low access
Counties immediately surrounding Philadelphia have higher percentage of population with low access to stores. In fact, relative to the rest of the state, Philadelphia has a lower percentage of individuals with low access to stores.
#Create new dataframe with variables of interest
LA2015PA <- Access %>%
select(FIPS, State, PCT_LACCESS_POP15) %>%
filter(State=="PA")
LA2015PA$GEO_ID <- paste0("0500000US", LA2015PA$FIPS)
LA2015PAmap <- inner_join(counties, LA2015PA, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorLA2015PAmapplot <- ggplot() +
geom_sf(data = LA2015PAmap) +
geom_sf(data = LA2015PAmap, aes(fill = PCT_LACCESS_POP15), lwd = 0) +
my_theme() +
ggtitle("County-Level Percent Population with Low Access to Stores in Pennsylvania, 2015") +
scale_fill_gradientn(name = "% of population with low access to stores", colours = myPalette(100))
print(LA2015PAmapplot)
Low Access and Race: Philadephia, PA versus Baltimore, MD
Surprisingly, although the overall food insecurity rate in Baltimore is lower than in Philadelphia, the percentage of individuals with low access to stores is much higher in Baltimore than in Philadelphia. This is likely due to the fact that Baltimore has far fewer grocery stores than Philadelphia, and instead has a larger proportion of convenience stores.
Philadelphia and Baltimore do not vary very much in terms of racial differences in access to stores.
#Create dataframe of interested variables
Access_race_PHLBAL <- Access %>%
select(FIPS, County, PCT_LACCESS_POP15, PCT_LACCESS_WHITE15,
PCT_LACCESS_BLACK15,
PCT_LACCESS_HISP15,
PCT_LACCESS_NHASIAN15,
PCT_LACCESS_NHNA15,
PCT_LACCESS_NHPI15,
PCT_LACCESS_MULTIR15) %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(Overall = PCT_LACCESS_POP15, White = PCT_LACCESS_WHITE15, Black = PCT_LACCESS_BLACK15, Hispanic = PCT_LACCESS_HISP15,
Asian = PCT_LACCESS_NHASIAN15,
American_Indian_Alaskan = PCT_LACCESS_NHNA15,
Hawaiian_Pacific_Islander = PCT_LACCESS_NHPI15,
Multiracial = PCT_LACCESS_MULTIR15) %>%
select(-FIPS)
ARPB <- data.frame(county = rep(c("Baltimore", "Philadelphia"), 8),
Race = rep(c("Overall", "White", "Black", "Hispanic", "Asian", "American_Indian_Alaskan", "Hawaiian_Pacific_Islander", "Multiracial"), each=2),
Percent = c(Access_race_PHLBAL$Overall, Access_race_PHLBAL$White, Access_race_PHLBAL$Black, Access_race_PHLBAL$Hispanic, Access_race_PHLBAL$Asian, Access_race_PHLBAL$American_Indian_Alaskan, Access_race_PHLBAL$Hawaiian_Pacific_Islander, Access_race_PHLBAL$Multiracial))
#Create a bar graph to visualize
ARPBplot <- ggplot(data = ARPB, aes(x=Race, y=Percent, fill = county)) + #Load data, specify variables
geom_bar(stat="identity", position = "dodge") + #Add a visual layer that is a barplot
scale_x_discrete(limits=c("Overall", "White", "Black", "Hispanic", "Asian", "American_Indian_Alaskan", "Hawaiian_Pacific_Islander", "Multiracial")) + #To manually adjust the order of the bars on the graph
labs(title="Relationship Between Race and Access to Stores in Philadelphia and Baltimore 2015") + labs(x="Race", y="Percent (%) Population with Low Access to Stores") +
theme(axis.text.x=element_text(angle=45,hjust=0.5,vjust=0.5))
ggsave(plot=ARPBplot, filename=paste("/Users/nawarnaseer/BMIN503_Final_Project/Plot_Images/", "ARPBplot.png", sep=""), width=7, height=5, units="in")Low Access and Age: Philadephia, PA versus Baltimore, MD
In Philadelphia, seniors have slightly better access than children to stores, whereas the opposite holds true for Baltimore.
#Create dataframe of interested variables
Access_age_PHLBAL <- Access %>%
select(FIPS, County, PCT_LACCESS_POP15,
PCT_LACCESS_CHILD15,
PCT_LACCESS_SENIORS15) %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(Overall = PCT_LACCESS_POP15, Children = PCT_LACCESS_CHILD15, Seniors = PCT_LACCESS_SENIORS15) %>%
select(-FIPS)
AAPB <- data.frame(county = rep(c("Baltimore", "Philadelphia"), each=3),
Age = rep(c("Overall", "Child", "Senior"), 2),
Percent = c(17.5, 3.98, 2.52, 1.21, 0.158, 0.221))#Create a bar graph to visualize
AAPBplot <- ggplot(data = AAPB, aes(x=Age, y=Percent, fill = county)) + #Load data, specify variables
geom_bar(stat="identity", position = "dodge") + #Add a visual layer that is a barplot
scale_x_discrete(limits=c("Overall", "Child", "Senior")) + #To manually adjust the order of the bars on the graph
labs(title="Relationship Between Age and Access to Stores in Philadelphia and Baltimore 2015") + labs(x="Age Group", y="Percent (%) Population with Low Access to Stores")
print(AAPBplot)
Financial access to food
We will first look at the overlap between individuals with low income and low access to stores. Next, we will look at price of a basic food item, milk, across the United States as an indirect measure of affordability of food, and thus financial access. The dataset also contains information about food assistance programs such as the Supplemental Nutrition Assistance Program (SNAP) and the Special Supplemental Nutrition Program for Women, Infants, and Children (WIC). We will use SNAP- and WIC- authorized stores as another indirect measure of financial access to food.
Low Access to Stores and Low Income
The map below displays the percent population with low access to stores that also had low income in 2015. This map tracks similarly to the map for low access, as does the map specifically for PA.
#Create new dataframe with variables of interest
LAI2015 <- Access %>%
select(FIPS, State, PCT_LACCESS_LOWI15)
LAI2015$GEO_ID <- paste0("0500000US", LAI2015$FIPS)
LAI2015map <- inner_join(counties, LAI2015, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vectorLAI2015mapplot <- ggplot() +
geom_sf(data = LAI2015map) +
geom_sf(data = LAI2015map, aes(fill = PCT_LACCESS_LOWI15), lwd = 0) +
my_theme() +
ggtitle("County-Level Percent Population with Low Income and Low Access to Stores in the US in 2015") +
scale_fill_gradientn(name = "Percent Population", colours = myPalette(100), limit=range(0, 100))
print(LAI2015mapplot)
PA-specific low access to stores and low income
LAI2015mapPAplot <- ggplot() +
geom_sf(data = LAI2015mapPA) +
geom_sf(data = LAI2015mapPA, aes(fill = PCT_LACCESS_LOWI15), lwd = 0) +
my_theme() +
ggtitle("County-Level Percent Population with Low Income and Low Access to Stores in Pennsylvania in 2015") +
scale_fill_gradientn(name = "Percent Population", colours = myPalette(100))
print(LAI2015mapPAplot)
Price of milk
The map below shows the price of low-fat milk as a ratio of the national average in 2010. Ratios greater than 1 indicate that the price is higher than the national average, while ratios lower than 1 indicate that the price is lower. This data is only available at the regional level. The price is higher in the southern and eastern parts of the US. The high price in the southern US may contribute to the high food insecurity in this region. The limitation of this dataset is that it does not include other commodities, and the price of low-fat milk is not sufficient to extrapolate prices of other goods.
#Create new dataframe with variables of interest
Price$GEO_ID <- paste0("0500000US", Price$FIPS)
Pricemap <- inner_join(counties, Price, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Price <- paste0(Pricemap$County, ",", " ", Pricemap$State, # paste0 to append tract name with other relevant text
"<br>Price of low-fat milk/national average, 2010: ", round(Pricemap$MILK_PRICE10, 4))
# Basic leaflet map
leaflet(Pricemap) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(MILK_PRICE10), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Price) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~MILK_PRICE10, # variable to be passed to palette function
title = 'Price of low-fat milk/national average, 2010', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()SNAP-authorized stores
#Create new dataframe with variables of interest
SNAP2016 <- Stores %>%
select(FIPS, State, County, SNAPSPTH16)
SNAP2016$GEO_ID <- paste0("0500000US", SNAP2016$FIPS)
SNAP2016map <- inner_join(counties, SNAP2016, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_SNAP2016 <- paste0(SNAP2016map$County, ",", " ", SNAP2016map$State, # paste0 to append tract name with other relevant text
"<br>SNAP-authorized stores/1,000 population, 2016: ", round(SNAP2016map$SNAPSPTH16, 1))
# Basic leaflet map
leaflet(SNAP2016map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(SNAPSPTH16), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_SNAP2016) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~SNAPSPTH16, # variable to be passed to palette function
title = 'SNAP-authorized stores/1,000 population, 2016', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific SNAP-authorized stores
In Pennsylvania, Philadelphia has the highest number of SNAP-authorized stores relative to the population.
#Create new dataframe with variables of interest
SNAP2016mapPA <- SNAP2016map %>%
filter(State=="PA")
SNAP2016mapPAplot <- ggplot() +
geom_sf(data = SNAP2016mapPA) +
geom_sf(data = SNAP2016mapPA, aes(fill = SNAPSPTH16), lwd = 0) +
my_theme() +
ggtitle("SNAP-authorized stores/1,000 population, 2016") +
scale_fill_gradientn(name = "SNAP-authorized stores/1,000 population, 2016", colours = myPalette(100))
print(SNAP2016mapPAplot)
WIC-authorized stores
There is a greater density of WIC-authorized stores in counties in the northern US.
#Create new dataframe with variables of interest
WIC2012 <- Stores %>%
select(FIPS, State, County, WICSPTH12)
WIC2012$GEO_ID <- paste0("0500000US", WIC2012$FIPS)
WIC2012map <- inner_join(counties, WIC2012, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_WIC2012 <- paste0(WIC2012map$County, ",", " ", WIC2012map$State, # paste0 to append tract name with other relevant text
"<br>WIC-authorized stores/1,000 population, 2012: ", round(WIC2012map$WICSPTH12, 1))
# Basic leaflet map
leaflet(WIC2012map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(WICSPTH12), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_WIC2012) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~WICSPTH12, # variable to be passed to palette function
title = 'WIC-authorized stores/1,000 population, 2012', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific WIC-authorized stores
In Pennsylvania, Philadelphia has the highest number of WIC-authorized stores relative to the population.
#Create new dataframe with variables of interest
WIC2012mapPA <- WIC2012map %>%
filter(State=="PA")
WIC2012mapPAplot <- ggplot() +
geom_sf(data = WIC2012mapPA) +
geom_sf(data = WIC2012mapPA, aes(fill = WICSPTH12), lwd = 0) +
my_theme() +
ggtitle("WIC-authorized stores/1,000 population, 2012") +
scale_fill_gradientn(name = "WIC-authorized stores/1,000 population, 2012", colours = myPalette(100))
print(WIC2012mapPAplot)
SNAP- and WIC-authorized Stores: Philadelphia, PA versus Baltimore, MD
Philadelphia has far higher numbers of SNAP- and WIC-authorized stores compared to Baltimore.
#Create dataframe of interested variables
SNAPWIC_PHLBAL <- Stores %>%
select(FIPS, County, SNAPSPTH12, WICSPTH12) %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(SNAP = SNAPSPTH12, WIC = WICSPTH12) %>%
select(-FIPS)
SWPB <- data.frame(county = rep(c("Baltimore", "Philadelphia"), 2),
Type = rep(c("SNAP", "WIC"), each=2),
Number = c(SNAPWIC_PHLBAL$SNAP, SNAPWIC_PHLBAL$WIC))
#Create a bar graph to visualize
SWPBplot <- ggplot(data = SWPB, aes(x=Type, y=Number, fill = county)) + #Load data, specify variables
geom_bar(stat="identity", position = "dodge") + #Add a visual layer that is a barplot
labs(title="Presence of SNAP/WIC-authorized Stores in Philadelphia and Baltimore 2012") + labs(x="Type of Assistance", y="Number of authorized stores per 1000 of population")
print(SWPBplot)
Utilization of food
This dataset does not contain any purchase or consumption data. However, we can use the number of SNAP/WIC participants and redemption of SNAP/WIC benefits as a measure of utilization of these benefits. Since both of these benefit programs have restrictions on the types of food items that can be purchased with them, and are generally limited to healthier, unprepared food items, we are using these as indirect measures for purchase and consumption of healthy food.
SNAP Participants
The map below shows the percentage of eligible individuals who participate in SNAP at the state level. Generally, there appears to be higher participation in the eastern United States. There does not appear to be a large difference between Maryland (90%) and Pennsylvania (90.4%).
#Create new dataframe with variables of interest
SNAPPar2016 <- Assistance %>%
select(FIPS, State, County, SNAP_PART_RATE13)
SNAPPar2016$GEO_ID <- paste0("0500000US", SNAPPar2016$FIPS)
SNAPPar2016map <- inner_join(counties, SNAPPar2016, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_SNAPPar2016 <- paste0(SNAPPar2016map$County, ",", " ", SNAPPar2016map$State, # paste0 to append tract name with other relevant text
"<br>SNAP Participants (% of Population), 2013: ", round(SNAPPar2016map$SNAP_PART_RATE13, 1))
# Basic leaflet map
leaflet(SNAPPar2016map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(SNAP_PART_RATE13), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_SNAPPar2016) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~SNAP_PART_RATE13, # variable to be passed to palette function
title = 'SNAP Participants (% of Eligible Population), 2013', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()SNAP Redemption
Overall, there is higher redemption of SNAP benefits in the western United States, the opposite of where there are higher number of participants. There is not much of a difference in redemption between Baltimore (approximately $250,000) and Philadelphia (approximately $360,000).
#Create new dataframe with variables of interest
SNAPRed2016 <- Assistance %>%
select(FIPS, State, County, REDEMP_SNAPS16)
SNAPRed2016$GEO_ID <- paste0("0500000US", SNAPRed2016$FIPS)
SNAPRed2016map <- inner_join(counties, SNAPRed2016, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_SNAPRed2016 <- paste0(SNAPRed2016map$County, ",", " ", SNAPRed2016map$State, # paste0 to append tract name with other relevant text
"<br>$ SNAP Redemption/SNAP-authorized stores, 2016: ", round(SNAPRed2016map$REDEMP_SNAPS16, 1))
# Basic leaflet map
leaflet(SNAPRed2016map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(REDEMP_SNAPS16), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_SNAPRed2016) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~REDEMP_SNAPS16, # variable to be passed to palette function
title = 'SNAP Redemption/SNAP-authorized stores, 2016', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific SNAP redemption
In Pennsylvania, Philadephia has high utilization of SNAP.
#Create new dataframe with variables of interest
SNAPRed2016mapPA <- SNAPRed2016map %>%
filter(State=="PA")SNAPRed2016mapPAplot <- ggplot() +
geom_sf(data = SNAPRed2016mapPA) +
geom_sf(data = SNAPRed2016mapPA, aes(fill = REDEMP_SNAPS16), lwd = 0) +
my_theme() +
ggtitle("SNAP Redemption, Pennsylvania, 2016") +
scale_fill_gradientn(name = "$ SNAP Redemption/SNAP-authorized stores, 2016", colours = myPalette(100))
print(SNAPRed2016mapPAplot)
WIC Partipants
The map below shows the percentage of individuals who participate in WIC at the state level. California and Texas stand out as having more participants than any other state. However, it should be noted that this data is not a percentage of the eligible population, but the population of the state as a whole. WIC is specifically for women, infants, and children, and there may simply be more eligible individuals in California and Texas relative to other states.
#Create new dataframe with variables of interest
WICPar2015 <- Assistance %>%
select(FIPS, State, County, PCT_WIC15)
WICPar2015$GEO_ID <- paste0("0500000US", WICPar2015$FIPS)
WICPar2015map <- inner_join(counties, WICPar2015, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_WICPar2015 <- paste0(WICPar2015map$County, ",", " ", WICPar2015map$State, # paste0 to append tract name with other relevant text
"<br>WIC Participants (% of Population), 2015: ", round(WICPar2015map$PCT_WIC15, 1))
# Basic leaflet map
leaflet(WICPar2015map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(PCT_WIC15), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_WICPar2015) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~PCT_WIC15, # variable to be passed to palette function
title = 'WIC Participants (% of Population), 2015', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()WIC Redemption
It is difficult to make any interpretations or conclusions due to the high number of missing values for this dataset.
#Create new dataframe with variables of interest
WICRed2012 <- Assistance %>%
select(FIPS, State, County, REDEMP_WICS12)
WICRed2012$GEO_ID <- paste0("0500000US", WICRed2012$FIPS)
WICRed2012map <- inner_join(counties, WICRed2012, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_WICRed2012 <- paste0(WICRed2012map$County, ",", " ", WICRed2012map$State, # paste0 to append tract name with other relevant text
"<br>$ WIC Redemption/WIC-authorized stores, 2012: ", round(WICRed2012map$REDEMP_WICS12, 1))
# Basic leaflet map
leaflet(WICRed2012map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(REDEMP_WICS12), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_WICRed2012) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~REDEMP_WICS12, # variable to be passed to palette function
title = 'WIC Redemption/WIC-authorized stores, 2012', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific WIC redemption
WIC utilization is not as high in Philadelphia relative to its immediate surrounding counties.
#Create new dataframe with variables of interest
WICRed2012mapPA <- WICRed2012map %>%
filter(State=="PA")WICRed2012mapPAplot <- ggplot() +
geom_sf(data = WICRed2012mapPA) +
geom_sf(data = WICRed2012mapPA, aes(fill = REDEMP_WICS12), lwd = 0) +
my_theme() +
ggtitle("WIC Redemption, Pennsylvania, 2012") +
scale_fill_gradientn(name = "$ WIC Redemption/WIC-authorized stores, 2012", colours = myPalette(100))
print(WICRed2012mapPAplot)
WIC Redemption: Philadelphia, PA versus Baltimore, MD
There is a significant difference in WIC redemption between Baltimore and Philadelphia. However, this could be due to differences in the percentage of eligible population in these two counties, which we do not have data on.
#Create dataframe of interested variables
WICRed2012_PB <- WICRed2012 %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(WIC = REDEMP_WICS12) %>%
select(-FIPS)#Create a bar graph to visualize
WICRed2012_PBplot <- ggplot(data = WICRed2012_PB, aes(x=County, y=WIC)) + #Load data, specify variables
geom_bar(stat="identity") + #Add a visual layer that is a barplot
labs(title="WIC Redemption in Philadelphia and Baltimore 2012") + labs(x="County", y="$ WIC Redemption / WIC-authorized Stores")
print(WICRed2012_PBplot)
Health outcomes
Food insecurity can have adverse effects on various health outcomes. In this section, we explore spatial distribution of obesity and diabetes to assess any overlaps with the distribution of food insecurity.
Diabetes
Diabetes is more prevalent in the southeastern part of the US, which does overlap with areas of high food insecurity.
#Create new dataframe with variables of interest
Diabetes2013 <- Health %>%
select(FIPS, State, County, PCT_DIABETES_ADULTS13)
Diabetes2013$GEO_ID <- paste0("0500000US", Diabetes2013$FIPS)
Diabetes2013map <- inner_join(counties, Diabetes2013, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Diabetes2013 <- paste0(Diabetes2013map$County, ",", " ", Diabetes2013map$State, # paste0 to append tract name with other relevant text
"<br>Adult Diabetes Rate (%), 2013: ", round(Diabetes2013map$PCT_DIABETES_ADULTS13, 1))
# Basic leaflet map
leaflet(Diabetes2013map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(PCT_DIABETES_ADULTS13), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Diabetes2013) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~PCT_DIABETES_ADULTS13, # variable to be passed to palette function
title = 'Adult Diabetes Rate, 2013', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific diabetes
In Pennsylvania, prevalence of diabetets is higher in Philadelphia relative to surrounding counties.
#Create new dataframe with variables of interest
Diabetes2013mapPA <- Diabetes2013map %>%
filter(State=="PA")Diabetes2013mapPAplot <- ggplot() +
geom_sf(data = Diabetes2013mapPA) +
geom_sf(data = Diabetes2013mapPA, aes(fill = PCT_DIABETES_ADULTS13), lwd = 0) +
my_theme() +
ggtitle("Adult Diabetes Rate, 2013") +
scale_fill_gradientn(name = "% of adults with diabetes, 2013", colours = myPalette(100))
print(Diabetes2013mapPAplot)
Obesity
Obesity is more prevalent in the eastern part of the US. It’s interesting to note that the same part of the United States with high rates of food insecurity have high rates of obesity.
#Create new dataframe with variables of interest
Obesity2013 <- Health %>%
select(FIPS, State, County, PCT_OBESE_ADULTS13)
Obesity2013$GEO_ID <- paste0("0500000US", Obesity2013$FIPS)
Obesity2013map <- inner_join(counties, Obesity2013, by="GEO_ID")## Warning: Column `GEO_ID` joining factor and character vector, coercing into character vector# Pop-up message
pu_Obesity2013 <- paste0(Obesity2013map$County, ",", " ", Obesity2013map$State, # paste0 to append tract name with other relevant text
"<br>Adult Obesity Rate (%), 2013: ", round(Obesity2013map$PCT_OBESE_ADULTS13, 1))
# Basic leaflet map
leaflet(Obesity2013map) %>%
addPolygons(stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun2(PCT_OBESE_ADULTS13), # addPolygons requires variable and a pass to a palette function
fillOpacity = 0.5, smoothFactor = 0.5, # increase opacity and resolution
popup = pu_Obesity2013) %>% # add a popup message
addTiles() %>% # default basemap
addLegend("bottomright", # location of legend
pal=pal_fun2, # palette function
values=~PCT_OBESE_ADULTS13, # variable to be passed to palette function
title = 'Adult Obesity Rate, 2013', # legend title
opacity = 1) %>% # legend opacity (1 = completely opaque)
addScaleBar()PA-specific obesity
Prevalence of obesity higher in Philadelphia relative to surrounding counties.
#Create new dataframe with variables of interest
Obesity2013mapPA <- Obesity2013map %>%
filter(State=="PA")Obesity2013mapPAplot <- ggplot() +
geom_sf(data = Obesity2013mapPA) +
geom_sf(data = Obesity2013mapPA, aes(fill = PCT_OBESE_ADULTS13), lwd = 0) +
my_theme() +
ggtitle("Adult Obesity Rate, 2013") +
scale_fill_gradientn(name = "% of adults with obesity, 2013", colours = myPalette(100))
print(Obesity2013mapPAplot)
Diabetes and Obesity: Philadelphia versus Baltimore
Overall rates of obesity and diabetes are both slightly higher in Philadelhia than Baltimore. Both counties have much higher rates of obesity than diabetes.
#Create dataframe of interested variables
Health_PHLBAL <- Health %>%
select(FIPS, County, PCT_DIABETES_ADULTS13, PCT_OBESE_ADULTS13) %>%
filter(FIPS %in% c("42101", "24005")) %>%
rename(Diabetes = PCT_DIABETES_ADULTS13, Obesity = PCT_OBESE_ADULTS13) %>%
select(-FIPS)
HPB <- data.frame(county = rep(c("Baltimore", "Philadelphia"), 2),
Outcome = rep(c("Diabetes", "Obesity"), each=2),
Percent = c(Health_PHLBAL$Diabetes, Health_PHLBAL$Obesity))#Create a bar graph to visualize
HPBplot <- ggplot(data = HPB, aes(x=Outcome, y=Percent, fill = county)) + #Load data, specify variables
geom_bar(stat="identity", position = "dodge") + #Add a visual layer that is a barplot
labs(title="Health Outcomes 2013") + labs(x="Health Outcome", y="Percent Rate")
print(HPBplot)
Conclusion
The maps generated in this report provide useful insight into spatial differences in factors influencing food insecurity throughout the United States. While this report focused on certain geographical areas (e.g. Philadelphia, PA and Baltimore, MD), the analyses generated in this report can be further studied and applied to other states and counties as well. The USDA Food Atlas also contains several other variables and indicators that were not incorporated into this report, but that would be useful to investigate in future studies. While the Food Atlas is very comprehensive in terms of variables, one limitation is that not all variables have data available at the county level. Additionally, not all variables have data collected from the same year. As such, comparisons of different variables across time was not possible. It would be interesting to incorporate data from other sources (e.g. NHANES, CDC datasets, etc.) for future analyses. Additionally, this report focused exclusively on the United States. Comparison of food insecurity in other developed countries may provide useful information on how best to design effective interventions.
Acknowledgements
I would like to thank the following individuals for their guidance and support throughout the generation of this report.
The following faculty/staff were consulted for guidance on designing the research question and selecting an appropriate dataset:
- Allison Karpyn, PhD
- Moriah Hall, MPH
- Ricardo Castillo-Neyra, PhD, DVM, MSPH
- Blanca Himes, PhD
The following individuals were instrumental in understanding statistics and troubleshooting errors with R code chunks:
- Erin Schnellinger, MS
- Emma Zheng, MS
- Jimson Mathew, BS, MA
References
- Coleman-Jensen, A., Rabbitt, M. P., Gregory, C. A., & Singh, A. (2019, September). Household Food Security in the United States in 2018. Retrieved November 6, 2019, from https://www.ers.usda.gov/publications/pub-details/?pubid=94848.
- Hunger Free America. Philadelphia Falling Behind in Freedom to Eat: Greater Philadelphia Hunger Report, 2018.
- Centers for Disease Control and Prevention, Leading Causes of Death, https://www.cdc.gov/nchs/fastats/leading-causes-of-death.htm
- Feeding America. Help Kids Facing Hunger This School Year, https://www.feedingamerica.org/hunger-blog/help-kids-facing-hunger-this